library(dplyr)
suppressPackageStartupMessages(library(GEOquery)) 
#gset=AnnoProbe::geoChina('GSE11318')
load("data/GSE11318_eSet.Rdata")
library(AnnoProbe)
library(limma)

set.seed(1024)

#基因表达矩阵


get_express_matrix <- function(gset){
  probes_expr <- gset[[1]] %>%  exprs()
  probes_expr <- probes_expr+1  %>% log2()
}

# 临床表型数据

get_pho_data <- function(gset){
  pdata <- pData(gset)
}

# 获取基因注释信息及探针

get_gpl_probe <- function(probes_expr,gset){
  probe2gene <- gset@annotation %>% idmap()
  genes_expr <- filterEM(probes_expr,probe2gene)
}

# make group list 为差异基因分析做准备
# 在此之前，进行抽样，使得数据平衡
# 在group里面抽样使得数据平衡
sample_balanced_inside <- function(clinical_beataml_BM,types,ratio = 1){
  long_df <- clinical_beataml_BM %>% filter(group == types[1]) 
  short_df <- clinical_beataml_BM %>% filter(group == types[2]) 
  
  #sam <- sample(1:nrow(short_df),nrow(long_df))
  #short_df <- short_df[sam,]
  if (nrow(long_df) >= nrow(short_df)){
    sam <- sample(1:nrow(long_df),nrow(short_df)*ratio)
    long_df <- long_df[sam,]
  }else{
    sam <- sample(1:nrow(short_df),nrow(long_df)*ratio)
    short_df <- short_df[sam,]
  }
  df <- rbind(long_df,short_df)
}

get_group_samples <- function(gset,types){
  pdata <- pData(gset) %>% select(characteristics_ch1.6,geo_accession) %>% 
    filter(characteristics_ch1.6 != "Clinical info: Final microarray diagnosis: Unclassified DLBCL") %>% 
    mutate(group = case_when(
      characteristics_ch1.6 == "Clinical info: Final microarray diagnosis: PMBL" ~ "PMBL",
      characteristics_ch1.6 == "Clinical info: Final microarray diagnosis: GCB DLBCL" ~ "GCBDLBCL",
      characteristics_ch1.6 == "Clinical info: Final microarray diagnosis: ABC DLBCL" ~ "ABCDLBCL"
    ) ) %>% 
    filter(group %in% types) %>% 
    sample_balanced_inside(types,ratio = 1) %>% 
    rename(sample_ID = geo_accession) %>% 
    select(-characteristics_ch1.6) %>% 
    relocate(sample_ID) %>% 
    arrange(group)
}

get_group_design <- function(group_sample){
  group_list=factor(group_sample$group)
  table(group_list)
  design=model.matrix(~factor(group_list))
  return(design)
}

get_expression_sort <- function(gset,group_sample){
  probes_expr <- exprs(gset)
  probes_expr <- probes_expr+1  %>% log2()
  
  probe2gene <- gset@annotation %>% idmap()
  genes_expr <- filterEM(probes_expr,probe2gene) %>% select(all_of(group_sample$sample_ID))
  
  return(genes_expr)

}
# 进行deg分析

perform_DEG_analysis <- function(genes_expr,design){
  genes_expr %>% 
    lmFit(design) %>% 
    eBayes() %>% 
    topTable(coef=2,n=Inf)
}

# DEG workflow

DEG_workflow <- function(gset,types){
  
  group_sample <- gset %>% get_group_samples(types = types)
  design <- group_sample %>% get_group_design()
  DEG <- gset %>% 
    get_expression_sort(group_sample) %>% 
    perform_DEG_analysis(design)
  return(DEG)
}

## visualization

deg_volcano_graph <- function(DEG,id){
  #DEG <- DEG[[ids]]
  need_deg=data.frame(symbols=rownames(DEG), logFC=DEG$logFC, p=DEG$P.Value)
  deg_volcano(need_deg,id)
}

deg_heatmap_graph <- function(DEG,genes_expr,group_list,id){
  #DEG <- DEG[[ids]]
  #genes_expr <- genes_expr[[ids]]
  #group_list <- group_list[[ids]]
  deg_heatmap(DEG,genes_expr,group_list,id)
}

# 正式计算

type_list <- list(
                  type1 = c("PMBL","ABCDLBCL"),
                  type2 = c("PMBL","GCBDLBCL"))
group_sample <- list()
genes_expr <- list()
DEG <- list()
gene_list <- list()
#vocalno_graphs <- list()
#hemaps <- list()
#check_diff_ccl5 <- list()
#check_diff_ccr5 <- list()
for (i in 1:2){
    DEG[[i]] <- DEG_workflow(gset[[1]],types = type_list[[i]])
    group_sample[[i]] <- gset[[1]] %>% get_group_samples(types = type_list[[i]])
    genes_expr[[i]] <-  gset[[1]] %>% 
      get_expression_sort(group_sample[[i]])
    gene_list[[i]] <- factor(group_sample[[i]]$group)
    #vocalno_graphs[[i]] <- deg_volcano_graph(DEG,i,2)
    #hemaps[[i]] <- deg_heatmap_graph(DEG,genes_expr,group_list,i,30)
    #check_diff_ccl5[[i]] <- check_diff_genes('CCL5',genes_expr[[i]],group_list[[i]])
    #check_diff_ccr5[[i]] <- check_diff_genes('CCR5',genes_expr[[i]],group_list[[i]])
    
}

# PMBL v.s. ABCDLBCL

ccl5_DEG_pa <- DEG[[1]]
groulist <- gene_list[[1]]
genesexp <- genes_expr[[1]] 
# 火山图
deg_volcano_graph(ccl5_DEG_pa,2)
# 热图
deg_heatmap_graph(ccl5_DEG_pa,genesexp,groulist,30)
# CCL5的差异检验
library(ggpubr)
data.frame(value=as.numeric(genesexp["CCL5",]),
              group=groulist) %>% 
ggboxplot( "group", "value",
          color = "group", palette =c("#00AFBB", "#E7B800"),
          add = "jitter", shape = "group")
ccl5_DEG_pa["CCL5",]
# CCR5的差异检验
data.frame(value=as.numeric(genesexp["CCR5",]),
           group=groulist) %>% 
  ggboxplot( "group", "value",
             color = "group", palette =c("#00AFBB", "#E7B800"),
             add = "jitter", shape = "group")
ccl5_DEG_pa["CCR5",]

# PMBL v.s. GCBDLBCL
ccl5_DEG_pa2 <- DEG[[2]]
groulist2 <- gene_list[[2]]
genesexp2 <- genes_expr[[2]] 
# 火山图
deg_volcano_graph(ccl5_DEG_pa2,2)
# 热图
deg_heatmap_graph(ccl5_DEG_pa2,genesexp2,groulist2,30)
# CCL5的差异检验
check_diff_genes('CCL5',genesexp2,groulist2)
ccl5_DEG_pa2["CCL5",]
# CCR5的差异检验
check_diff_genes('CCR5',genesexp2,groulist2)
ccl5_DEG_pa2["CCR5",]
